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Abstract 


Research  on  combustion  is  being  conducted  at 
Lewis  Research  Center  to  provide  improved  analyti¬ 
cal  models  of  the  complex  flew  and  chemical  reac¬ 
tion  processes  which  occur  in  the  combustor  of  gas 
turbine  engines  ana  other  aeropropulsion  systems. 
The  objective  of  the  research  is  to  obtain  a  bet¬ 
ter  understanding  of  the  various  physical  proces¬ 
ses  that  occur  in  the  gas  turbine  combustor  in 
order  to  develop  models  and  numerical  codes  which 
can  accurately  describe  these  processes.  Activi¬ 
ties  include  in-house  research  projects,  univer¬ 
sity  grants,  and  industry  contracts  and  are 
classified  under  the  subject  areas  of  advanced 
numerics,  fuel  sprays,  fluid  mixing,  and  radiation- 
chemistry.  Results  are  highlighted  from  several 
projects. 


Introduction 

As  we  enter  into  the  age  of  the  computer  and 
the  laser,  new  opportunities  are  emerging  to  make 
significant  advances  in  the  understanding  of  the 
physics  of  combustion.  These  advances  are  needed 
in  the  gas  turbine  engine  industry.  Development 
costs  and  maintenance  costs  are  prevalent  factors 
when  considering  the  next  generation  of  gas  tur¬ 
bine  engines.  The  costs  associated  with  the 
design  and  development  of  engine  components, 
including  the  combustor,  are  alarmingly  high  using 
today's  design  methodology.  Also,  the  ever- 
increasing  operating  pressures  and  temperatures  of 
new  engines  for  increased  performance  and  better 
fuel  economy  are  having  adverse  effects  on  dura¬ 
bility  of  hot  section  components.  This  results  in 
high  maintenance  costs.  Clearly  what  is  needed  is 
an  improved  design  methodology  which  can  allow  the 
design  engineer  to  converge  to  the  optimum  design 
in  a  much  shorter  development  cycle. 

At  NASA  Lewis  Research  Center,  combustion 
research  is  integrated  into  activity  in  internal 
fluid  mechanics.  The  objectives  of  the  research 
are  to  advance  the  understanding  of  flow  physics, 
heat  transfer  and  combustion  processes  which  are 
fundamental  to  aeropropulsion,  and  to  translate 
this  knowledge  into  models  and  numerical  codes  of 
aerothermodynamic  phenomena.  The  overall  goal  is 
to  bring  internal  computational  fluid  mechanics  to 
a  state  of  practical  application  to  propulsion 
systems.  These  models  and  numerical  codes  would 
then  be  available  to  the  industry  to  incorporate 
into  their  own  engine/component  design  systems. 

The  approach  at  Lewis  is  to  establish  an 
integrated  computational-experimental  research 
program.  The  activity  consists  of  research  on 


numerical  methods,  well-defined  experiments  for 
code  and  model  verification,  and  the  demonstration 
of  computational  codes  for  propulsion  system 
components.  In  the  area  of  combustion,  research 
has  been  underway  for  some  time  in  these  three 
areas. ’*3  with  a  recent  reorganization  in  the 
Aeronautics  Program  at  Lewis,  this  research  activ¬ 
ity  has  been  consolidated  with  similar  research  in 
other  components  such  as  compressors,  turbines,  and 
transition  ducts  and  is  continuing  in  the  internal 
fluid  mechanics  research  program.  Although  the 
amount  of  work  currently  being  supported  in  com¬ 
bustion  has  been  decreased  in  the  new  organiza¬ 
tion,  the  prime  research  needs  in  combustion  are 
recognized  and  the  program  is  focusing  on  them. 


Figure  1  is  a  cut-away  view  of  a  ficticious 
combustor  which  illustrates  the  complex  fluid 
mechanics  and  combustion  features.  The  flows  are 
highly  three-dimensional  with  turbulence  levels  in 
many  cases  comparable  in  magnitude  to  the  bulk 
velocity.  Liquid  fuels  are  injected  as  a  spray 
which  then  undergoes  vaporization  and  mixing. 
Chemical  reaction  occurs  which  causes  changes  in 
density  and  fluid  mechanics  properties  and  can 
cause  the  formation  of  a  solid  phase  (soot)  with 
its  attend-.'t  high  radiatior  heat  transfer.  The 
understanding  of  these  physical  processes  is 
needed  before  accurate  numerical  codes  can  be 
built  and  used  as  a  predictive  tool  in  the  design 
process.  In  addition  the  numerical  methods  for 
three-dimensional  flows  need  improvements  in 
accuracy  and  efficiency  in  order  to  properly  simu¬ 
late  the  features  of  these  flows.  This  then  is 
the  framework  of  the  current  combustion  research 
program.  As  shown  in  Fig.  ?,  the  activities  are 
focused  in  four  areas:  advanced  numerics,  fuel 
sprays,  fluid  mixing,  and  radiation-chemistry. 

Both  experiments  and  computational  research  are 
being  conducted,  with  the  long  range  objective  of 
providing  the  numerical  codes  that  can  be  used 
with  confidence  as  a  predictive  tool  in  the  indus¬ 
try's  design  system.  In  this  report  a  research 
project  from  each  of  these  areas  will  now  be  high¬ 
lighted  as  examples  of  the  kinds  of  work  currently 
being  supported  at  Lewis. 


Advanced  Numerics:  Error  Reduction  Program 


A  major  restriction  to  the  development  of  a 
computer  based  design  methodology  is  the  accuracy 
of  the  numerical  methods  used  in  combustor  flow 
codes.  The  upwind  differencing  currently  used  in 
these  flow  codes,  introduces  an  appreciable 
error  in  the  calculated  result.  This  error  (or 
numerical  diffusion)  is  frequently  of  such  a  large 
magnitude  that  it  obscures  the  turbulence  model 
used  in  the  calculation. 


To  alleviate  this  problem,  NASA  has  conducted 
a  program  to  identify  and  incorporate  an  improved 
accuracy  differencing  scheme  into  a  combustor  flow 
code.  Under  a  portion  of  this  program  a  variety 
of  differencing  schemes  were  examined  in  several 
test  calculations.  The  schemes  examined  included 
QUICK  (Quadratic  Upstream  Interpolation)  and  SUD 
(Skewed  Upwind  Differencing) .  QUICK  differencing 
was  developed  by  Leonard.®  This  scheme  improves 
the  accuracy  of  convective  differencing  by  per¬ 
forming  an  upwind  biased  quadratic  interpolation. 
This  scheme  is  second  order  accurate  and  can  pro¬ 
duce  nonphysical  oscillations  in  the  solution. 

SUD,  developed  by  Raithby,®  attains  high  accu¬ 
racy  by  differencing  in  an  upwind  manner  along  the 
flow  streamlines.  While  maintaining  the  same  for¬ 
mal  accuracy  as  upwind  differencing  the  truncation 
error  in  SUD  is  smaller  in  magnitude.  As  with 
QUICK,  SUD  can  produce  nonphysical  oscillations  in 
the  solution,  therefore,  a  scheme  to  "bound"  SUD 
was  also  examined.  This  scheme  employs  the  con¬ 
cept  of  flux-blending/  wherein  a  bounded  flux 
determined  from  upwind  differencing  is  blended 
with  the  unbounded  -  but  more  accurate  -  SUD  flux. 
The  main  factor  is  to  blend  as  little  of  the  lesser 
accurate  scheme  while  still  maintaining  a  properly 
"bounded"  solution.  This  procedure  called  8SUDS, 
starts  from  an  initial,  totally  skew  differenced 
estimate  and  blends  an  upwind  flux  if  the  solution 
is  out  of  the  range  of  neighboring  values.  If  the 
solution  is  in  range,  i.e.,  bounded,  then  no  blend¬ 
ing  is  performed. 

An  illustration  of  the  accuracy  of  the 
upwind,  QUICK,  and  SUD  schemes  is  seen  in  Fig.  3. 
The  figure  displays  the  results  of  a  single  point 
scalar  transport  calculation  made  for  various  flow 
angles.  All  schemes  agree  with  an  exact  solution 
(no  error)  at  a  flow  angle  of  zero,  but  departing 
from  this  each  scheme  displays  some  degree  of 
error  relating  to  numerical  diffusion.  The  error 
displayed  by  upwind  differencing  increases  with 
flow  angle  up  to  a  maximum  at  45*.  QUICK  displays 
a  similar  behavior  but  the  overall  level  of  error 
is  much  less.  SUD  displays  an  error  maximum 
around  15*  but  tends  to  zero  at  angles  approaching 
45*.  Both  QUICK  and  SUD  display  a  much  higher 
level  of  accuracy  than  upwind. 

A  scalar  transport  calculation  is  useful  for 
generally  examining  some  aspects  of  differencing 
scheme  performance,  but  a  more  complete  test  is  a 
laminar  flow  calculation.  The  results  of  a  series 
of  laminar  flow  calculations  from  Ref.  8  are  dis¬ 
played  in  Fig.  4.  In  this  figure  axial  velocity 
profiles  at  a  distance  of  one-half  a  duct  height 
from  the  inlet  are  shown  for  two  different  compu¬ 
tational  meshes.  In  these  calculations  an  indica¬ 
tor  of  accuracy  is  the  steepness  of  the  velocity 
profile.  Steep  velocity  profiles  are  exhibited  by 
QUICK  and  BSUDS  with  the  upwind  profiles  exhibit¬ 
ing  a  high  degree  of  numerical  diffusion.  On  the 
coarse  mesh,  BSUDS  appears  to  be  more  accurate 
than  QUICK,  while  on  a  fine  mesh  there  is  not  much 
of  a  distinction  between  the  two  schemes. 

An  important  aspect  of  improved  accuracy  is 
the  effect  that  the  differencing  scheme  has  on  the 
rate  of  convergence.  The  computational  times  to 
converge  the  system  of  governing  equations  in  the 
previous  laminar  calculations  is  shown  in  Table  1. 

I  he  convergence  times  are  ratioed  to  the  upwind 
convergence  times  to  clearly  illustrate  the  compu¬ 
tational  penalty  paid  to  attain  improved  accuracy. 


Generally,  the  improved  accuracy  schemes  required 
from  3  to  15  times  longer  to  reach  a  converged 
solution.  To  a  degree  one  can  hope  that  this 
computational  penalty  can  be  offset  with  the  use 
of  coarse  meshes  to  achieve  the  same  overall  level 
of  accuracy.  In  reality  it  appears  that  a  rela¬ 
tively  fine  mesh  is  needed  even  with  the  high 
accuracy  schemes.  In  any  case  the  need  for 
improved  solution  algorithms  with  these  more  accu¬ 
rate  differencing  schemes  is  strongly  demonstrated. 

Fuel  Sprays:  Study  of  Dilute  Turbulent 
Particle-Laden  Jets 

A  research  topic  of  current  interest  in  com¬ 
putational  fluid  mechanics  is  spray  processes  in 
combustion  chambers.  However,  due  to  the  compli¬ 
cations  involved  in  the  spray  processes,  e.g., 
hydrodynamics  of  spray  atomizatior,  polydisperse 
drop  diameter  distributions,  drop  shattering  and 
coalescence,  and  evaporation  and  combustion  of 
drops,  spray  processes  in  combustion  chambers  are 
excessively  complex  for  systematic  model  develop¬ 
ment  and  evaluation.  In  contrast,  injection  and 
mixing  of  monodisperse  solid  particles  is  a  much 
simpler  process,  and  facilitates  the  modeling  of 
the  underlining  physics  involved  in  practical  mul¬ 
tiphase  flows. 

The  objective  of  this  work  was  to  evaluate 
models  of  dilute  turbulent  particle-laden  jets. 
Dilute  solid-particle-laden  jets  in  a  still  envi¬ 
ronment  were  considered.  These  flows  are  rela¬ 
tively  simple  and  only  involve  interphase  momentum 
exchange.  Modeling  efforts  concentrated  on  effects 
of  turbulent  fluctuations  on  momentum  transfer 
between  the  phases,  as  well  as  the  dispersion  of 
the  particles  by  turbulent  fluctuations.  Three 
typical  models  of  the  flow  processes  were  consid¬ 
ered:  (1)  a  locally  homogeneous  flow  (LHF)  model, 
where  slip  between  the  phases  was  neglected;  (2)  a 
deterministic  separated  flow  (DSF)  model,  where 
slip  was  considered  but  effects  of  turbulence  on 
particle  motion  were  ignored;  and  (3)  a  stochastic 
separated  flow  (SSF)  model  where  effects  of  both 
interphase  slip  and  turbulence  on  particle  motion 
were  considered  using  random-sampling  techniques. 

The  LHF  model  represents  the  simplest  treat¬ 
ment  of  a  multiphase  flow.  In  this  model,  the 
particles  and  the  continuous  phase  are  assumed  to 
have  equal  rates  of  turbulent  diffusion.  This 
implies  that  interphase  transport  rates  are  infi¬ 
nitely  fast,  so  that  both  phases  have  the  same 
velocity  at  each  point  in  the  flow.  The  LHF 
approximation  is  valid  only  for  flows  containing 
small  particles,  where  characteristic  response 
times  of  particles  are  small  in  comparison  to 
characteristic  times  of  turbulent  fluctuations. 

F aeth’  has  indicated  that  LHF  models  only  yield 
accurate  predictions  for  particle  sizes  smaller 
than  most  practical  applications.  Despite  the 
limitations,  LHF  models  have  some  important  advan¬ 
tages.  They  require  minimal  information  concern¬ 
ing  initial  conditions  of  particle  properties, 
which  are  often  difficult  to  obtain  for  practical 
flows.  The  theoretical  model  of  the  flow  is  equiv¬ 
alent  to  that  of  a  single-phase  flow  and  effects  of 
multiple  phases  only  appear  in  the  representation 
of  state  relationships,  bypassing  the  difficulties 
involved  with  modeling  interactions  between  the 
phases. 
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The  DSF  model  considers  slip  between  the 
phases  but  neglects  turbulent  particle  dispersion. 
In  this  model  particles  follow  deterministic  tra¬ 
jectories  and  only  interact  with  mean  properties 
of  the  continuous  phase.  The  computation  involves 
dividing  the  particulate  phase  into  representative 
samples  whose  motion  and  transport  are  tracked 
through  the  flow  field  using  a  Lagrangian  formula¬ 
tion.  An  Eulerian  formulation  is  employed  to 
solve  the  governing  equations  for  the  continuous 
phase.  The  effect  of  interphase  transport  is  con¬ 
sidered  by  inserting  appropriate  source  terms  in 
the  governing  equations  for  the  continuous  phase. 
The  source  terms  are  determined  using  the  "Particle 
Source  in  Cell"  technique.  Neglecting  turbulent 
dispersion  is  appropriate  for  flows  containing 
large  particles,  where  particle  response  times  are 
large  in  comparison  to  characteristic  turbulent 
fluctuation  times. 

Most  practical  particle-laden  flows  and 
sprays,  however,  exhibit  properties  between  the 
two  limits  represented  by  LHF  and  OSF  models  and 
therefore  require  consideration  of  turbulent  par¬ 
ticle  dispersion.  The  stochastic  separated-f low 
(SSF)  approach,  first  proposed  by  Gosman  and 
Ioannides*'*  and  subsequently  developed  by  Faeth 
and  coworkers, was  adopted  to  study  the 
effects  of  turbulence  on  particle  dispersion.  The 
SSF  model  involves  finding  trajectories  of  statis¬ 
tically  significant  samples  of  individual  parti¬ 
cles  as  they  move  away  from  the  injector  and 
encounter  a  random  distribution  of  turbulent 
eddies.  Particle  trajectory  computations  are 
similar  to  the  DSF  model,  except  that  instanta¬ 
neous  eddy  properties  replace  mean-gas  properties. 
Eddy  properties  are  obtained  by  constructing  the 
probability  density  function  (PDF)  of  velocity, 
assuming  isotropic  turbulence,  and  random  sampling 
to  find  each  velocity  component. 

The  three  models  of  particle-laden  jets  have 
been  systematically  evaluated  and  some  typical 
results  are  shown  here.  Comparisons  of  model 
predictions  with  experimental  data  are  shown  in 
Figs.  5  and  6,  taken  from  Ref.  13.  Figure  5 
presents  radial  profiles  of  mean  and  fluctuating 
properties  for  the  gas  phase,  and  Fig.  6  presents 
radial  profiles  for  the  particle  phase.  Since 
this  jet  has  a  relatively  low  particle  loading, 
gas  phase  properties  are  relatively  independent  of 
model  type  and  both  the  LFHF  and  SSF  models  are  in 
good  agreement  with  the  measurements.  For  the 
particle  properties  SSF  predictions  agree  reason¬ 
ably  well  with  the  data.  Predictions  of  the  LHF 
model  for  particle  properties  are  not  at  all 
satisfactory,  since  neglect  of  slip  overestimates 
particle  fluctuations  levels  and  rates  of  spread. 
The  DSF  model  also  yields  poor  results  because 
neglecting  particle  dispersion  by  turbulence 
causes  the  rate  of  spread  of  the  particles  to  be 
substantially  underestimated.  Similar  results 
were  found  throughout  the  data  base  used  in  the 
study. 

For  the  two  separated-f low  models  used  in 
this  study  (DSF  and  SSF)  the  model  predictions  are 
relatively  insensitive  to  the  specification  of 
gas-phase  initial  conditions.  The  specification 
of  initial  particle  properties,  however,  exerts 
much  more  pronounced  effects  on  predictions. 
Accurate  measurement  of  particle  initial  condi¬ 
tions,  e.g.,  particle  size,  mean  and  fluctuating 
velocities  and  particle  mass  flux,  is  essential  to 


obtain  meaningful  predictions  of  separated-f low 
models,  especially  when  particle  relaxation  time 
is  large  compared  to  particle  residence  time  in 
the  flow  field. 

Practical  particle-laden  flows  and  sprays 
usually  involve  regions  where  the  dilute  flow 
approximation  is  not  appropriate,  such  as  near  the 
injector.  Phenomena  to  be  considered  for  future 
model  extension  in  application  involving  dense 
particulate  flows  include:  the  finite  volume 
fraction  of  the  dispersed  phase,  effects  of  nearby 
particles  on  particle  transport  properties,  and 
effects  of  particle  collision  and  coalescence. 
There  is  also  a  pressing  need  to  overcome  the 
experimental  difficulties  encountered  in  the  dense 
flow  region  so  that  further  progress  can  be  made 
toward  developing  reliable  models  for  this  region. 

Fluid  Mixing:  Direct  Numerical  Simulations  of 
Chemically-Reacting  Mixing  Layers 

The  advent  of  high  speed  vector  processing 
computers  has  made  possible  the  study  of  turbu¬ 
lence  and  turbulence-chemistry  interactions  via 
computational  "experiments."  One  promising 
approach  that  has  been  supported  at  NASA  Lewis*4 
is  Direct  Numerical  Simulation  (DNS).  Direct 
Numerical  Simulations  solve  the  complete  time- 
dependent,  Navier-Stokes  equations  at  a  range  of 
flow  Reynolds  number  wherein  all  the  scales  of 
turbulence  can  be  fully  resolved  on  a  computa¬ 
tional  mesh.  This  typically  requires  the  calcula¬ 
tion  of  flows  where  the  Reynolds  number  based  on 
the  Taylor  macroscale  is  of  the  order  of  100.  The 
key  assumption  involved  is  that  the  large  scales 
of  turbulence  which  can  be  represented  on  the  com¬ 
putational  mesh  will  not  significantly  vary  with 
Reynolds  number.  While  this  may  not  be  true  for 
all  the  aspects  of  turbulent  flow  which  one  might 
desire  to  study,  a  previous  study  of  incompres¬ 
sible  turbulent  flows*^  has  indicated  that  a 
great  deal  of  insight  into  the  structure  of  turbu¬ 
lent  flows  can  be  obtained.  One  of  the  main  goals 
is  that  this  information  will  be  effective  in 
developing  turbulence-closure  models. 

Incompressible  turbulence  models  typically 
require  information  on  fluctuating  pressure- 
velocity  correlations  that  are  difficult  to 
measure  experimentally.  Reacting  flows  further 
complicate  the  problem  with  the  need  for  an 
instantaneous  measure  of  density  and  the  number  of 
double  and  triple  correlations  rapidly  becomes 
overwhelming.  Indeed,  even  the  form  of  averaging 
(Fauvre  or  Reynolds)  to  be  used  in  constructing 
the  governing  equations  is  not  straight-foward.*^ 
The  promise  of  DNS  is  that  it  can  provide  complete 
flow  field  values  that  can  be  processed  to  provide 
any  necessary  correlation  or  test  various  closure 
models. 

The  above  states  the  promise  of  DNS,  but 
there  are  a  number  of  "realities"  which  should 
also  be  mentioned.  The  first  has  already  been 
noted:  simulations  must  be  made  at  low  Reynolds 
number.  This  is  primarily  a  numerical  resolution 
limitation.  Even  employing  highly  accurate 
psuedo-spectral  methods*7  and  large  computa¬ 
tional  meshes  (typically  64  by  64  by  64  for  the 
results  reported  herein)  the  range  of  turbulent 
scales  that  can  be  represented  is  limited.  A 
second  limitation  is  that  the  numerical  methods 
used  restrict  the  calculation  to  relatively  simple 
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geometries.  The  numerical  calculations  discussed 
here  are  of  a  time-evolving  shear  layer  as  opposed 
to  the  experimentally  measured  spacial  mixing 
layer.  The  time-evolving  mixing  layer  is  not  a 
direct  analog  of  the  spacial  mixing  layer,  and 
some  physical  phenomenon  cannot  be  represented, 
tor  example,  the  flow  streamlines  in  the  time- 
•rvolving  case  are  not  aligned  along  the  flow 
direction  as  they  are  in  the  spacial  case.  A 
final  limitation  relates  to  the  statistical  nature 
of  DNS  results.  One  simulation  must  typically  run 
for  several  CPU  hours  to  simulate  a  small  fraction 
of  physical  time  (on  the  order  of  a  second).  Sta¬ 
tistics  sensitive  to  scale  variation  may  not  be 
aoequately  sampled.  It  is  easily  seen  that  more 
small  scales  of  turbulence  can  be  represented  on  a 
t  mite  computational  domain  than  larger  scales  and 
tins  can  bias  the  statistical  results.  These  fac¬ 
tors  require  care  in  the  sampling  of  simulation 
results  to  obtain  statistically  stationary 
1 nformation. 

Despite  the  current  limitations  inherent  to 
INS  a  number  of  turbulent  flow  phenomenon  can  be 
well  represented.  An  example  is  the  experimental 
versus  computational  flow  visualization  results 
shown  in  Fig.  7.  The  experimental  flow  visualiza¬ 
tion  was  provided  by  Koochesf ahani1®  and  dis¬ 
plays  a  flow  structure  resulting  from  the  merging 
of  two  discrete  vortices.  This  flow  structure  is 
strikingly  similiar  to  the  results  of  a  two- 
dimensional  numerical  simulation.  The  experiment 
was  conducted  at  a  Reynolds  number  and  with  fluids 
differing  from  the  computation,  yet  the  comparison 
is  encouraging. 

Another  favorable  aspect  of  DNS  results  is 
oispl’ayed  in  Fig.  8.  Similarity  theory  predicts 
that  the  mean  velocity  half-width  of  the  shear 
layer  should  grow  linearly.  As  seen  in  Fig.  8, 
the  time  evolution  of  the  velocity  half-width 
grows  linearly  with  time  in  the  numerical  simula¬ 
tions.  If  a  transformation  relating  space  and 
time  is  performed1'  the  time  evolving  results 
ure  found  to  be  in  good  agreement  with  similarity 
: henry.  This  agreement  is  favorable  for  a  large 
range  of  initial  conditions  as  shown  on  Fig.  8. 

A  final  comparison  between  simulated  results 
ana  experimental  data  is  shown  in  Fig.  9.  The 
experimental  data,  Mungal,1^  is  of  a  mildly  exo- 
Li.ermic  chemical  reaction.  The  similarity  plot  of 
integrated  product  level  from  a  numerical  simula¬ 
tion  agrees  fairly  well  with  the  shape  and  peak  of 
i >e  experimental  data.  The  shape  of  the  product 
profile  is  somewhat  sensitive  to  the  initial  con¬ 
ditions  used  in  the  simulation  and  current  research 
is  being  conducted  to  assess  whether  alternate 
initial  conditions  might  clear  up  the  small  dis¬ 
crepancy  displayed. 

The  overall  agreement  between  DNS  results  and 
experimental  information  is  strongly  encouraging. 

1  he  highlights  of  the  current  research  is  merely 
illustrated  here  and  Ref.  14  is  recommended  for 
further  details. 

From  this  array  of  favorable  initial  results 
it  appears  that  this  approach  can  provide  insights 
into  the  physics  of  turbulent  chemistry  which  will 
one  day  lead  to  improved  closure  models  usable  by 
combustor  designers. 


Radiation,  Chemistry:  Fast  Algorithms  for 
uombustion  Kinetics  tabulations 

The  ordinary  differential  equations  (ODE's) 
describing  complex  chemical  reactions  are  charac¬ 
terized  by  widely  different  time  constants. 

Although  the  ODE's  are  stable,  standard  numerical 
techniques  such  as  the  popular  explicit  Runge-Kutta 
and  Adams  methods  are  prohibitively  expensive  to 
use  because  of  the  severe  steplength  restrictions 
imposed  by  the  requirements  for  numerical  stabil¬ 
ity.^  Such  systems  of  ODE's  are  commonly 
referred  to  as  “stiff''  systems.21 

The  aim  of  the  present  investigation  was  to 
examine  currently  available  integration  techniques 
for  solving  chemical  kinetic  rate  equations.  The 
motivation  behind  this  work  is  the  increasing 
interest  in  (1)  modeling  the  reaction  mechanisms 
describing  the  consumption  of  fuels  and  pollutant 
formation  and  destruction,  and  (2)  multidimensional 
modeling  of  reactive  flows,  which  includes  the 
equations  of  fluid  motion.  The  former  results  in 
the  need  to  integrate  large  systems  of  nonlinear 
ODE's.  The  latter  results  in  the  need  to  inte¬ 
grate  the  ODE's  at  several  thousand  grid  points. 

To  make  such  calculations  practicable,  it  is 
necessary  to  have  a  very  fast  homogeneous  batch 
chemistry  integrator. 

The  techniques  examined  include  two  general- 
purpose  solvers,  EPISODE22  and  LSODE,  ®  developed 
for  an  arbitrary  system  of  ODE's,  and  the  special¬ 
ized  techniques  CHEMEQ,24  CREK1D,25  and  GCKP84,26 
developed  specifically  for  chemical  kinetics  appli¬ 
cations.  Testing  of  these  codes  was  accomplished 
by  application  to  two  practical  combustion  kinetics 
problems.  Both  problems  described  adiabatic,  homo¬ 
geneous,  gas-phase,  transient,  batch  combustion 
reactions  at  constant  pressure.  Test  problem  1 
described  the  ignition  and  subsequent  combustion 
of  a  mixture  of  33  percent  carbon  monoxide  and 
67  percent  hydrogen  with  100  percent  theoretical 
air.  It  consisted  of  12  reactions  among  11  spe¬ 
cies.  Test  problem  2,  involving  30  reactions  among 
15  species,  described  the  ignition  and  subsequent 
combustion  of  a  stoichiometric  hydrogen-air  mix¬ 
ture.  Both  problems  were  integrated  over  a  time 
period  of  1  msec,  which  encompassed  all  three  com¬ 
bustion  regimes:  induction,  heat  release,  and 
equilibration. 

In  using  LSODE,  EPISODE,  and  CHEMEQ  the  tem¬ 
perature  was  calculated  by  employing  two  different 
methods.  In  method  A  the  temperature  was  evalu¬ 
ated  via  an  iterative  solution  of  the  algebraic 
enthalpy  conservation  equation.  In  method  B  the 
temperature  was  computed  by  solving  its  ODE.  Both 
CREK1D  and  GCKP84  were  written  explicitly  for  non- 
isothermal  chemical  reactions  and  therefore  include 
calculation  procedures  for  the  temperature. 

Figures  10  and  11  present  the  variation  of 
the  local  error  tolerance  with  the  computer  ( i . e . , 
CPU)  time  required  to  solve  test  problems  1  and  2, 
respectively.  In  these  figures  CPU  is  the  com¬ 
puter  time  (in  seconds)  required  on  the  NASA  Lewis 
Research  Center's  IBM  370/3033  computer.  All 
results  were  obtained  using  single-precision 
accuracy,  except  GCKP84  which  was  in  double- 
precision.  The  error  control  performed  is  differ¬ 
ent  for  the  different  codes.  Consequently,  the 
local  error  tolerance,  EPS,  does  not  have  the  same 
meaning  for  all  the  codes.  For  LSODE,  GCKP84, 


CREKlD,  ana  CHEMEQ,  EPS  is  the  local  relative 
error  tolerance/^  For  EPISODE,  however,  EPS  is 
a  mixed  relative  and  absolute  error  tolerance  - 
relative  for  species  with  initially  nonzero  mole 
fractions  ( i • e • ,  reactants)  and  for  the  tempera¬ 
ture,  and  absolute  for  species  with  initially  zero 
mole  fractions  (i.e.,  intermediate  species  and 
products) . 

Figures  10  and  11  show  that  LSODE  is  the 
fastest  code  currently  available  for  solving  com¬ 
bustion  kinetic  rate  equations.  Because  GCKP84 
has  been  designed  for  a  variety  of  chemical  kinet¬ 
ics  problems  it  requires  more  work  per  step  than 
LSODE  and  is  hence  slower.  CREKlD  and,  especially 
for  test  problem  2,  EPISODE  are  attractive  alter¬ 
natives  to  LSODE.  However,  there  are  many  prob¬ 
lems  associated  with  using  EPISODE  for  solving 
combustion  kinetics  rate  equations  -  details  are 
given  in  Ref.  20.  Figures  10  and  11  show  also 
that  the  use  of  temperature  method  A  does  not 
result  in  significant  inefficiency.  On  the  con¬ 
trary  this  method  can  be  significantly  faster  than 
evaluating  the  temperature  by  integrating  its  ODE 
(Fig.  11). 

For  multidimensional  chemically  reacting  flow 
calculations  computational  speed  is  of  primary 
concern  and  moderate  accuracy  is  adequate.  How¬ 
ever,  in  developing  and  validating  reaction  mecha¬ 
nisms,  accuracy  is  of  critical  importance.  An 
accuracy  comparison  of  the  techniques  examined  in 
this  study  was  made  as  follows.  For  each  test 
problem  a  mean  integrated  root  mean  square  error 
( E rms)  was  estimated.  Erms  is  approximately 
the  average  true  error  incurred  by  a  technique  in 
solving  the  complete  problem.  Figures  12  and  13 
present  the  variation  of  Er(ns  with  the  user- 
supplied  value  for  the  local  error  tolerance, 

EPS.  Note  the  significant  discrepancies  between 
values  specified  for  EPS  and  the  errors  actually 
obtained.  For  both  problems  temperature  method  A 
is  as  accurate  as  temperature  method  B  -  in  many 
cases,  it  is  significantly  more  accurate.  LSODE 
is  the  most  accurate  code  for  solving  chemical 
kinetic  rate  equations.  However,  GCKP84,  CREKlD, 
ana  ChEMEQ-A  compare  favorably  with  LSODE  for 
small  values  of  EPS.  EPISODE  is  significantly 
less  accurate  than  the  other  codes  because  the 
error  control  used  in  it  is  inappropriate  for 
chemical  kinetics  applications.^® 
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TABLE  I.  -  RATIO  OF  CONVERGENCE  TIMES 
FOR  VARIOUS  DIFFERENCING  SCHEMES 
WITH  UPWIND  CONVERGENCE  TIMES 
USED  AS  THE  STANDARD 
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•  FULLY  3-DIMENSIONAL  FLOW  •  CHEMICAL  REACTION/HEAT  RELEASE 

•  HIGH  TURBULENCE  LEVELS  *  2  PHASE  W,TH  VAPORIZATION 

Figure  1.  -  Cutaway  illustration  of  a  portion  of  a  full  annular 
combustor  system,  with  major  features  highlighted. 


RESEARCH  LONG  RANGE  APPLICATION 

ELEMENTS  RESEARCH  OBJECTIVES  BY  INDUSTRY 


Figure  2.  -  Combustion  research  has  long  range  objectives  to  produce 
predictive  numerical  codes  that  have  practical  application  in  engine 
manufacturer's  combustion  design  system. 


ERROR,  PERCEN1 


Figure  3.  -  Results  of  a  scalar  transport  test 
calculation  for  various  differencing  schemes. 
Error  from  an  exact  solution  is  displayed 
versus  flow  angle. 


(b)  Fine  grid  (58x38)  calculations. 

Figure  4  -  Laminar  flow  test  calculations  comparing 
upwind,  QUICK,  and  BSUDS.  Inlet  flow  angle  of  25° 
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(a)  Computational  results. 


(b)  Experimental  flow  visualization. 


Figure  7.  -  Flow  structure  comparison  between  direct  numerical 
simulation  and  experimental  flow  visualization. 


Figure  8.  -  Mean  velocity  half-width,  Zm  (■  Zj^- 
versus  time  -  three-dimensional  simulations. 
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Figure  9.  -  Similarity  plot  of  integrated  product 
level  versus  radial  location.  3-D  simulation. 
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Figure  10.  -  Variation  of  the  CPU  time(s)  with  error 
tolerance,  EPS,  for  test  problem  1.  All  runs  on 
IBM  370/3033. 
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Figure  11.  -  Variation  of  the  CPU  time(s)  with  error 
tolerance,  EPS,  for  test  problem  2.  All  runs  on 
IBM  370/3033. 
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